knitr::opts_chunk$set(echo = TRUE)
library(cytosignal)
library(ggplot2)

Load the data and plot False Positive Rates from real and simulated replates.

diff.data = readRDS("FP-diffusion-real_simulated.rds")
cont.data = readRDS("FP-contact-real_simulated.rds")

# We only have 5 real replicatesm, but 10 simulations
# Create the boxplot
p.diff <- ggplot(diff.data, aes(x = analysis, y = fp_rate, fill = analysis)) +
  # geom_boxplot() +
  geom_violin() +
  # add jitter points
  # geom_jitter(width = 0.01, height = 0.01, alpha = 0.3, size=0.1, color = "black") +
  scale_x_discrete(drop = FALSE)+
  # facet_wrap(~Metric, scales = "fixed") +
  ylim(0, 1) +
  cowplot::theme_cowplot(16)+
  theme(legend.position = "none")+
  # theme(axis.text.x = element_text(angle = 30, hjust = 0.9, vjust = 0.8))+
  xlab("Median of ligand/receptor gene expression in each bin") + ylab("Diffusion-dependent")

# same for the other lists
p.cont <- ggplot(cont.data, aes(x = analysis, y = fp_rate, fill = analysis)) +
  # geom_boxplot() +
  geom_violin() +
  # add jitter points
  # geom_jitter(width = 0.01, height = 0.01, alpha = 0.3, size=0.1, color = "black") +
  scale_x_discrete(drop = FALSE)+
  # facet_wrap(~Metric, scales = "fixed") +
  ylim(0, 1) +
  cowplot::theme_cowplot(16)+
  theme(legend.position = "none")+
  # theme(axis.text.x = element_text(angle = 30, hjust = 0.9, vjust = 0.8))+
  xlab("Median of ligand/receptor gene expression in each bin") + ylab("Contact-dependent")

all.data = rbind(diff.data, cont.data)

p.all <- ggplot(all.data, aes(x = analysis, y = fp_rate, fill = analysis)) +
  geom_violin() +
  scale_x_discrete(drop = FALSE)+
  ylim(0, 1) +
  cowplot::theme_cowplot(16)+
  theme(legend.position = "none")+
  xlab("Median of ligand/receptor gene expression in each bin") + ylab("Combined")

# Save the plots for reference
cowplot::plot_grid(p.diff, p.cont, p.all, nrow=3)

Load the data and plot False Negative Rates from two different perturbation simulations.

diff.data = readRDS("FN-diffusion-perturb.rds")
cont.data = readRDS("FN-contact-perturb.rds")

# We only have 5 real replicatesm, but 10 simulations
# Create the boxplot
p.diff <- ggplot(diff.data, aes(x = analysis, y = fn_rate, fill = analysis)) +
  geom_boxplot() +
  # geom_violin() +
  # add jitter points
  # geom_jitter(width = 0.01, height = 0.01, alpha = 0.3, size=0.1, color = "black") +
  scale_x_discrete(drop = FALSE)+
  # facet_wrap(~Metric, scales = "fixed") +
  ylim(0, 1) +
  cowplot::theme_cowplot(16)+
  theme(legend.position = "none")+
  # theme(axis.text.x = element_text(angle = 30, hjust = 0.9, vjust = 0.8))+
  xlab("Analysis") + ylab("Diffusion-dependent")

# same for the other lists
p.cont <- ggplot(cont.data, aes(x = analysis, y = fn_rate, fill = analysis)) +
  geom_boxplot() +
  # geom_violin() +
  # add jitter points
  # geom_jitter(width = 0.01, height = 0.01, alpha = 0.3, size=0.1, color = "black") +
  scale_x_discrete(drop = FALSE)+
  # facet_wrap(~Metric, scales = "fixed") +
  ylim(0, 1) +
  cowplot::theme_cowplot(16)+
  theme(legend.position = "none")+
  # theme(axis.text.x = element_text(angle = 30, hjust = 0.9, vjust = 0.8))+
  xlab("Analysis") + ylab("Contact-dependent")

all.data = rbind(diff.data, cont.data)

p.all <- ggplot(all.data, aes(x = analysis, y = fn_rate, fill = analysis)) +
  # geom_violin() +
  geom_boxplot() +
  scale_x_discrete(drop = FALSE)+
  ylim(0, 1) +
  cowplot::theme_cowplot(16)+
  theme(legend.position = "none")+
  xlab("Analysis") + ylab("Combined")

# Save the plots for reference
cowplot::plot_grid(p.diff, p.cont, p.all, nrow=3)

Load the data and plot the spatial heatmap colored by real data, regular simulation, and two perturbed simulations.

knitr::opts_chunk$set(echo = TRUE)
library(SingleCellExperiment)
library(ggplot2)
load("sce_reg_simu.RData")

plotGenes <- function(genes, sce, title){

  loc = colData(sce)[,c("x","y")]
  expre = lapply(genes, function(x){
    curr = as.matrix(counts(sce)[x,])
    curr = log1p(curr)
    # return(scales::rescale(curr))
    return(curr)
  })
  long = do.call(rbind, expre)
  long = as.data.frame(long)
  colnames(long) <- "Expression"
  long$gene = do.call(c, lapply(genes, function(x){rep(x,dim(expre[[1]])[1])}))
  long$x = rep(loc[,1], length(genes))
  long$y = rep(loc[,2], length(genes))

  p0 = as_tibble(long, rownames = "Cell") %>% ggplot(aes(x = x, y = y, color = Expression)) +
    geom_point(size = 0.1, stroke = 0.1)+facet_grid(~gene)+
    scale_colour_gradientn(colors = viridis_pal(option = "magma")(10), limits=c(0, 1)) +
    coord_fixed(ratio = 1) + theme(axis.text.x = element_text(angle = 45)) +
    cowplot::theme_cowplot(12)+
    ggtitle(title)

  return(p0)
}

de.genes = "BMP4"

de.real = plotGenes(de.genes, sce, "Real DE genes")
de.reg = plotGenes(de.genes, reg_sce, "Reg Simulated DE genes")
de.simu.p30k = plotGenes(de.genes, simu_sce.p30k, "Pertube 30K cells")
de.simu.p1k <- plotGenes(de.genes, simu_sce.p1k, "Pertube 1K cells")

# remove all x and y axis labels and titles
de.real <- de.real + theme(legend.position = "none",axis.text.x = element_blank(), axis.text.y = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank())
print(de.real)

de.reg <- de.reg + theme(legend.position = "none",axis.text.x = element_blank(), axis.text.y = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank())
de.reg

de.simu.p30k <- de.simu.p30k + theme(legend.position = "none",axis.text.x = element_blank(), axis.text.y = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank())
de.simu.p30k

de.simu.p1k <- de.simu.p1k + theme(legend.position = "none",axis.text.x = element_blank(), axis.text.y = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank())
de.simu.p1k